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Abstract 

The master equation for a charged harmonic oscillator coupled to an electromagnetic reservoir 
is investigated up to fourth-order in the interaction strength by using Krylov averaging method. 
The interaction is in the velocity-coupling form and includes a diamagnetic term. Exact analytical 
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environment. It is found that, generally, the third and the fourth order contributions have opposite 
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I. INTRODUCTION 


Research activity in dissipation and decoherence of quantum systems has accelerated 
in the last two decades with developments in quantum information technologies, progress 
in detection and control of ultra-cold atom and ions, atom lasers and photonic band-gap 
materials as well as the research on the role of decoherence in quantum-to-classical transi¬ 
tion [1-5]. Theoretical efforts to investigate dissipation and decoherence phenomena due to 
interaction between a quantum system and its environment use either exact or perturbative 
treatments. While exact approaches can be used to treat a restricted class of problems for 
the whole interaction parameter range, the perturbative treatments can be applied to almost 
all problems, but only for a restricted range of parameters which is, generally, the weak cou¬ 
pling regime. A plethora of different perturbation based formalisms have been developed to 
treat such system-bath problems over the last six decades [3-5]. The most prominent among 
them are Nakajima-Zwanzig projector operator method, path-integral based influence func¬ 
tional method [6], self-consistent hybrid schemes, Monte Carlo methods and hierarchical 
approaches [7, 8]. Most of the perturbative approaches express the equation of motion of 
the reduced density matrix of the system as a quantum master equation (QME) which 
contains second-order interaction terms. 

One way to increase the validity range of the perturbative methods might be increasing 
the order of master equation. Actually, it is shown by Fleming and Cummings [9] that an 
order-2n accuracy in the full-time solutions of master equations requires an order-(2n -|- 2) 
master equation. A number of groups investigated the effects of the fourth and higher order 
terms in detail [8-22]. In particular, Reichman et. al. have derived analytical expressions 
for dynamics of a harmonic oscillator in contact with various types of baths in the low tem¬ 
perature limit [15]. Also, a physically sound fourth-order correction to Redheld equation has 
been worked out by Laird, Budimir and Skinner [11]. Jang, Cao and Silbey [19] have derived 
a fourth-order QME in both time local and nonlocal forms for a general system Hamiltonian 
and used it to study motion of a particle in a continuous potential held and the dynamics 
of a two-level system coupled to a bath. Fourth-order corrections were found to cause a 
potential renormalization for the hrst problem while they introduced additional coherence 
for the later. In a series of papers, Brener et. al. [16-18] have examined, in detail, the effects 
of including the fourth and higher order terms in time convolutionless master equation for 
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the damped Jaynes-Cummings model (JCM), atom laser and the damped harmonic oscilla¬ 
tor. It was found that at weak and moderate coupling, the fourth order QME was in good 
agreement with the exact solution for the JCM while for the intermediate coupling regime 
of atom lasers, one needs to include non-Markovian as well as the fourth- and sixth-order 
terms to the master equation to adequately describe the dynamics of the system. Singh 
and Brumer investigated the validity of the of the second-order Markovian Redfield theory 
for reorganization energy and decay rates of photon autocorrelation functions in dimer sys¬ 
tems [22]. Liu et. al. [8] used hierarchical equation of motion method for the spin-boson 
problem and showed that fourth-order corrections are important for the intermediate cou¬ 
pling regime. Mavros and Voorhis considered fourth-order corrections to the memory kernel 
of the generalized master equation of the spin-boson problem and found that a numerically- 
exact solution was possible when the system-bath coupling is weak [23]. Interestingly, it was 
found that the non-perturbative treatment of certain system-bath problems is equivalent 
to the second-order perturbative master equations and including higher than second order 
terms would be detrimental for such systems [20, 21]. 

Dissipative dynamics of a charged particle in an electromagnetic held, whether a co¬ 
herent laser held or an incoherent blackbody radiation, have been examined by several 
groups [14, 24-32] because of its relevance to damped transport in solid-state systems. Bao 
and Bai [28, 29] have considered the the motion of a single electron atom interacting with an 
electromagnetic held in the dipole coupling approximation and obtained a ballistic dihusion 
for the particle in the position space. Kalandarov et. al. [30] have shown that dissipation of 
collective energy of a charged oscillator linearly coupled to a heat bath could be controlled 
by an external held. Pachon and Brumer [31] investigated the ehect of blackbody radiation 
on the coherence of the dynamics of light-harvesting molecules in photosynthesis and showed 
that the blackbody radiation would not enhance the coherence of the system. 

Ford, Lewis, and O’Connell derived a master equation for a charged oscillator weakly 
coupled to an electromagnetic held by using Krylov averaging method which can be improved 
by including higher-order interaction terms in a systematical way [14]. The model has several 
interesting features, such as use of velocity-coupling for the system-environment interaction 
which is not as widely studied as the position-coupling and keeping the diamagnetic term 
of the interaction Hamiltonian which leads to the nonzero third order contributions to the 
master equation which are absent for most of similar studies. In the current work, our 
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aim is two-fold; we aim for exact analytical expressions for all coefficients of the derived 
master equation which are obtained in the limit of infinite environmental cutoff frequency 
in Ref. [14], Our findings indicate that the exact treatment is crucial for the accurate high 
frequency behavior of the diffusion coefficients. Second purpose of the present study is 
to derive analytical expressions for the fourth-order corrections to the master equation of 
the charged oscillator in an electromagnetic field and delineate the convergence range of 
interaction parameters. 

The structure of this paper is as follows. We will introduce the model of the charged 
harmonic oscillator in contact with a blackbody radiation field in Section II where we pro¬ 
vide a short review of the derivation of the master equation by Krylov averaging method. 
Analytical expressions for the second, the third and the fourth order corrections are, also, 
presented in Section II and discussed in Section III. Section IV concludes the paper with 
a short review of the main findings. We provide the details of the derivation of the fourth 
order terms in the Appendix. 


II. MODEL AND METHOD 


We consider a charged harmonic oscillator in a blackbody radiation field [14, 33, 34]. The 
total Hamiltonian of the closed system formed by the oscillator and field can be written as: 

H = Hs + Hr + eHi, (1) 


with 


Hs = 
Hr = 


eHj = 



\kx\ 

— + vf), 

rrij 2 J 

2M’ 


( 2 ) 

( 3 ) 

( 4 ) 


where Hs, Hr, and eHj describe the charged harmonic oscillator as the system, blackbody 
radiation field as the reservoir and the interaction between the velocity of the oscillator and 
the vector potential A of the field, respectively. K and M are the force constant and the 
mass of the charged oscillator, while rrij, pj, qj and Uj describe the environmental oscillators 
that stand for the blackbody field, p and x for the oscillator and pj, qj of the reservoir obey 
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the usual commutation rules as [x,p\ = ih, [qi,Pj] = ih5ij and [x,pj] = [p,qj]=0. The vector 
potential of the held A can be expressed as: 


^ = ( 5 ) 

j 

The interaction between the oscillator and the held is assumed to be in velocity coupling 
type and contains both p.A and a diamagnetic A"^ terms. It is well known that the Hamil¬ 
tonian in Eq. (1), can be transformed to a renormalized coordinate-coupling form by using 
Power-Zienau transformation [35]. The diamagnetic interaction A^ is generally neglected 
[36, 37], because its inclusion makes calculation of path-integral based master equation dif- 
hcult. But its inclusion is shown to be important for the correct derivation of the partition 
function of a charged oscillator in blackbody held [34, 38] and we keep it here. 

The inhuence of the environment on the system is determined by the spectral density 
J (ca) of the bath which is a measure of the number of bath oscillators at frequency u and the 
coupling constant of the system-bath interaction. For the present case, we use the Drude 
corrected Ohmic spectral density: 


where Te = 2e^/(3mc^) = 6.24 10“^^ s is the 2/3 of the time it takes light to traverse the 
classical radius of electron. O is the cutoh frequency of the bath modes which sets bath 
correlation time as tb = 0“^. There is an upper limit of the form O < based on the 
causality arguments [39]. 

In the interaction picture, the unitary dynamics of the closed system formed by the 
oscillator and the radiation field is described by the von Neumann equation for the density 
matrix p as: 


where the interaction picture transformation is taken with respect to Hs and Hr as 


(7) 


eHi it) 

A{t) 


ijHs+Hjip i{Hg+Hn)t 

= e ft eii/e ft 
Pit) A it) , A^jt) 
m 2m 

= e ft H e ft 


(8a) 


= '^^{rrij Uj qj cos {uj t) + Pj sin {uj t)}, 

j 

p{t) = e ft pe ft 


(8b) 

(8c) 
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We assume that the bath and the oscillator are uncorrelated initially, that is 

P (0) = Ps (0) ® Pi?, (9) 

where ps = Tri?{p} is the reduced density matrix of the oscillator and pr is the density 
matrix of the reservoir which is assumed to be thermal at temperature T. Tri? indicates 
partial trace over the bath degrees of freedom. 

There have been may attempts to solve Eq. (7) by using various different approximations, 
such as projector operator [40, 41], path integral influence functional formulations [42], 
When the system-environment interaction is weak, Krylov averaging method can also be 
used to obtain a hierarchy of equations for the dynamics of the oscillator [14], For the self¬ 
containment of the present paper, we give a short overview of the method and the details of 
the derivation of the fourth order terms in the Appendix section. 

For ease of comparison with Ref. [14], we will rearrange the commutators on the right- 
hand side of Eq. (A7) as 

^ ^ ^ ^0 As [p (t ), [x (t ), Ps]] 

+ A 4 ([p(t), {pit) -imuoxit)}ps] + [psipit) +imuoxit)},pit)])) , ( 10 ) 


by dehning 

Ai = —Im (Ti), A 2 = — (Re (Ti) + Im (T 2 )), 

Uq Uq 

A3 = —Re(T2), A4 = —Im(T2). (11) 

Uq Uq 

At the second order, we need to evaluate and dehned in equations (A 8 a) and 
(A 8 b). Both of them can be calculated exactly by using the rational expansion of the 
hyperbolic cotangent 

coth ilSuo) = ^ 02 2 (^ 2 ) 

p Uq “ p^ ujn + 'n 

as follows: 


Ti = fuo (^7 coth if3uo) -iuo- 

Ti'^ = /.o7(/-0, 


(13b) 


where the mass renormalization factor 6m/m and the decay constant 7 are dehned in Ref. [14] 


= n Te and 7 = Wq Te , 
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and 


/ = 


1 2 


Re ipo 


1-i 


■ l3u)o 


n 


- ^0 




IT 


(15) 


where ipo [a)] is the digamma function and Re stands for the real part of its argument. One 

(o\ fo\ 

should note that T{ and T 2 of equations (A 8 a) and (A 8 b) correspond to C and S integrals 
of Ref. [14], respectively. They were evaluated at the O —)■ cx) limit in Ref. [14], which might 
cause some inconsistencies [43]. Using the mass renormalization and damping constants as 
dehned in Ref. [14], one can see that the main modihcation from the exact integration is 
scaling of these quantities by the electron structure factor = Vl?/{u^ + O^). scaling 

factor here is important because although mass renormalization is obtained by inhnite O 
limit, O needs to be hnite for the dehnition of 5m to be meaningful and f^og provides a 
cutoff for the calculated quantities at high oscillator frequency Uq, also. / term in Eq. (13b), 
which does not contribute to the master equation at the rotating wave approximation level, 
is also evaluated exactly here in terms of digamma functions and depends on the ratio of the 
temperature to the cutoff frequency of the blackbody held as well as the oscillator frequency. 

Using the exact integrals of and given by equations (13a) and (13b) in Eq. (11), 
we get the second order coefficients of master Eq. (10) 


m 0)0 




/I _ f J— 

^4 — Jujo 


(16) 


Wo Wo 

where N{I3uq) is the thermal occupation at inverse temperature /9 and oscillator frequency 


Uo. 


The third order contributions which are due, partly, to the diamagnetic term in the 
interaction Hamiltonian will be obtained by evaluating complex T 31 and T 32 and real T 33 and 
T 34 integrals dehned in the Appendix including equations (A9a)-(A9d). These integrals can, 
also, be evaluated exactly by using the rational expansion of coth(/9a;o) function Eq. (12). 
The results are 


T31 = 
T 32 = 
T 33 = 
T34 = 


I — coth {13 Uo) - — I 

m Uq 


- 7/1 


IwlQ 


7 , / ^ N 5m ^ 

— coth {f3 Uo) H- I 

Uo m 


, -o n2 

+ ^27 —A , 


Uq 


m 


- 7/1 


5m 


LJO 


1 


m \f3uo 


+ coth {/3 Uo) j , 


-7/010 — ^0 [/ 5 ,^,c^o], 

n m 



(17a) 

(17b) 

(17c) 

(17d) 
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where I is defined in Eq. (15), iho [P, i^o] = '^o ~ [l ~ } and ipo{x) is the 

digamma function. The dominant term in all four integrals involve product of decay rate 
and the mass renormalization and their squares which are modulated by the cutoff factor 
fojQ- The real parts of all four are temperature dependent while the imaginary parts of T 31 
and T 32 depend on only the mass renormalization and the interaction strength terms. One 
should note that, if the cutoff frequency of the blackbody field is chosen as the causally 
allowed highest value as O = then ^ = 1. Noting that ^ ~ 10“^ for the optical 
frequencies, one can drop the Tq [/S, O, wq] term of T 31 and coth {/3uo) term of T 32 . 

The third order contribution to the coefficients of Eq. (10) can be obtained by plugging 
the results displayed in equations (17a)-(17d) into the definitions in Eq. (11), which results 









Uo 




[l + 2iV(a;o)]-—/ 

Uq 


7 Sm 

+ jLJO 

ujq m 



1 

f3uo 



1_ 


^[l + 2iV(/3a;o)] 

cuo 


2 /. 


2 7 6m 
m 


+ l 1)1. 

m \ TT UoJ 1 


(18a) 

(18b) 

(18c) 

(18d) 


The most important difference between equations (18a)-(18d) and the third order coefficients 
in equation (7.16) of Ref. [14] is the scaling which provides frequency cutoff for each one 
of those coefficients. We, also, provide exact expressions for and integrals of Ref. [14] 
which were unevaluated there. 

The fourth order contributions are tedious but straightforward to obtain by performing 
the integrals defined in Eq. (A17) and are given in equations (Al 8 a)-(Al 8 e) of the Appendix. 



Here, we collect the fourth-order contributions to the coefficients of the Eq. (10): 




A ?=-/’ 


+ TT — 3 — ) ) [1 + 2 N {/3 (Uo). 



(19a) 


(19b) 


4 ( 4 ) - SL 

- Jujo,, 
Uo 


^6m 7 

m n 


(19c) 
(19d) 

where is the polygamma function of order i and Ini{'^i(/9, wq)} = lull'll(1— i/9 cuo/tt)}. 
All four coefficients have an overall scaling and are function of various combinations of 
interaction strength 7 /a;o, mass renormalization 6m/m and 7 /fl. and are inde¬ 

pendent of the temperature while and A^^^ have complicated temperature dependence 
through thermal occupation factor as well as various polygamma functions. 


III. DISCUSSION 

We will discuss the convergence of the expansion and the relative importance of various 
terms in Eq. (10) by a slight change of the master equation to 

^ ^ ~2^ ^ ^ 

+imuo\ \p, {x, ps}]+m uq D^p [p, [x, ps]]} , (20) 

where A = is mass renormalization, A = A^*^ is decay constant, D^x = — 

A 4 ^) is normal diffusion and Dxp = X^jAg^ is anomalous diffusion coefficients. To discuss 
the relevance of various terms in Eq. (20), we will state the adjoint master equation that 
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governs the dynamics of the moments of position and momentnm operators of the oscillator 
as [3] 

I (O), = {[Hs, 0|), + ([p^ O] >, - D„ (Ip. IP. 0]]>, 

+imuo\{{x, [p,0]})^ - muoD^p{[p, [x,0]])^}, ( 21 ) 


where D^x term causes diffusion of the variance of the oscillator position operator, Dxp plays 
the same role for x-p, A term is a change in the mass and A is decay constant of the position. 
It is obvious that while diffusion constants depend on the temperature, drift terms, mass 
renormalization and the damping constants, are temperature independent. 

By rearranging equations (16), (18a) and (19a) the renormalization constant A can be 
expressed in terms of Te, Uq and hi in a more suggestive form as 


A 


— |l + 



l) + fuQ (l — 3a;o) + oJq (l 



( 22 ) 


where uq = uq/VL and VL = VLXe are the dimensionless oscillator and the bath cutoff frequency, 
respectively and = 1/(1 + Wg) is the bath cutoff function. Note that the inverse of Uq 
can be considered as a resonance factor, its a measure of the effectiveness of the interaction 
between the charged oscillator and the environment, a large or small value of Aq indicates 
that the oscillator is detuned from the peak of the environmental spectral distribution. The 
successive terms in the nested parenthesis of the expression for A come from the second, the 
third and fourth order corrections, respectively. These terms are in the form of product of 
mass renormalization and cutoff factor modulated with the ratio of the oscillator frequency to 
the bath cutoff. The ratio of the third and the fourth order to the second order contribution 
can be expressed as (/(Ag — 1)/(Aq + 1) and 12^(1 — 3Aq)/(Aq + 1)^ which approach 12 and 0 
at the high frequency limit uq S> 12, respectively. In the low frequency limit Aq 1, these 
ratios tend to —12 and 12^. 

Similar to A, we can rewrite A in a simpler form as 


— —fcdo ^ 0 1^ 1 ^ + /lio^ + fcjQ (12(3 — Ag) + 2A| 


(23) 


The perturbative terms enter into the expression of A as powers of /lipH. A comparison of 
equations ( 22 ) and (23) shows that the leading terms of A and A are fcj^ 12 and fcjQ Aq 12, 
respectively. So, when the oscillator is in resonance with the peak frequency of the environ¬ 
ment, the magnitude of renormalization and the damping constants would be comparable 
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while off-resonance A is expected to be small compared to A for low frequencies. At high 
frequencies, —)■ 0 and A and A tend to zero. 



UJq j 

FIG. 1. Magnitude (a) and ratio of the third (b) and fourth order contributions (c) to the second 
order term for A and 7 as function of the oscillator frequency wq and the cutoff frequency of the 
blackbody environment fl. Both axis scales and the contour separations are logarithmic and the 
dashed lines indicate negative contours. 

The magnitude of A and A along with the ratio of the third and the fourth order contri¬ 
butions to their leading term, as function of Gre and ojq/VL, are presented in Fig. l(a)-(c) 
and Fig. l(d)-(f), respectively. Both the renormalization and the decay constant are nega¬ 
tive for the considered frequency range. A*^^^/A and changes sign at wq = H while 

A^^^A < 0 and A^^^^A > 0 for the whole range of Wq/H and Vt considered here. These ratios 
can be considered as a proxy for the convergence of the perturbative expansion of the master 
equation. It is clear that and contributions partially cancel each other and only 
the second order term B^‘^'> might be sufficient to account for the dynamics of the oscillator 
for the most part of the considered parameter range. Only very near the upper limit of the 


11 

































cutoff frequency (f2 ^ the whole perturbative treatment might be divergent because 
the magnitude as well as the contributions approach —1 and ±1, respectively. It is probable 
that including even higher order terms in the master equation will not lead to a converged 
A and A in this parameter regime. One should note that the interaction strength in the 
current model is y/tuo = ojQTf, and the non-convergence issue is independent of its value for 
high cutoff O. We display the uq dependent A and A at three different Ore values in Fig. (2) 
which indicates that while and contributions are signihcant at Ore = 0.5, they 
decrease substantially when Ore is reduced to 0.1 and become insignihcant at Ore = 10“^. 
It can be deduced from Fig. 1(b), (c), (e) and (f) that the relative contributions of the 
third and the fourth order terms are less than 1% for Ore < 0.01 and including only the 
second order terms in master equation might be sufficient for an accurate description of the 
dynamics of the charged oscillator. 


xIO"' xIO"' x10'^ x10'^ xio'^ xio'^ 



FIG. 2. The second (straight), the third (dot-dashed) and the fourth (dotted) order contributions 
to mass renormalization (left axis) and damping rate (right axis) as function of the dimensionless 
oscillator frequency at three different environmental cutoff frequencies (a) fire = 0.5, b) fire = 0.1 
and c) fire = 0.001). 

Both the normal and the anomalous diffusion coefficients are temperature dependent. We 
will hrst discuss the high temperature limit of the higher order contributions to D^x and Dxp- 
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Defining z = I3VL^ z corresponds to the high temperature regime and Taylor expanding 
hyperbolic cotangent and polygamma functions around zero and taking the hrst term of 
the expansion as, coth(a;o;2) —)■ 1/{Cjqz), cschcao;^ —)■ 1/CjqZ, '^i(1 ± iCjqz/ti) —)■ 7r^/6, 
'ilji{z/n) —)■ !z^ and 'ilj2{z/7r) —27r^jz^ we get 


£)( 2 ) 

XX 

£^(3) 

XX 

££(4) 

XX 


h - 

° 2 : ’ 


02 
f3 _ 
J ion . 


Z L 


3S) + (i^(2 - f!) 


(24) 


and 


££( 2 ) 

xp 

££(3) 

xp 

££(4) 

xp 


fcjQ Wo 




-2/|.c3„l^(2 + i,„-), 


4o ^0 


D 2 


L 


6D +6^2(3+ ca2)(l + D) 


(25) 


which indicate that all contributions to D^x and Dxp are linear in temperature at the high 
temperature limit. Both the second and the fourth order contributions to both diffusion 
coefficients are positive while the third order ones are always negative. 

One should note that from a comparison of equations (24)-(25) that the terms contribut¬ 
ing to Dxp contain an uq factor which is a measure of interaction strength. So, anomalous 
diffusion coefficient is expected to be much smaller than the normal diffusion coefficient in 
the weak coupling limit. In this limit (cuo —)■ 0), /(wq) = 1 and the convergence of the 
coefficients depends solely on the cutoff frequency of the environment as can be seen from 
equations (24)-(25). 

For the low temperature limit 3> 1, we can use the approximations coth (; 2 ) —)■ 1, 
csch (; 2 ) —)■ 0, Re{'^o(l — * Wo ^/tt) —ipo{z/n)} —)■ log(a;o), z\v[i[%Iji{ 1 — ioj^z/Ti)] —)■ tt/cDo, 
z‘^'ip 2 {z/TT) and zipi^z/n) ^ tt to express Dxx and Dxp as 




2 + Wo-Wo log (wo) , 

71 

I 2 “ ~ 2 ^ ^ + Wq + 2 log (cuo)] I , 


( 26 ) 
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and 


=-Uo^o^^og (wo), 
71 




UJQ 


(Do H—(2 + (Dq) log ((Do) 


TT 


—~ f^o "I 512 — 1 + (Dq( 2 + 12) + (1 + 612)(Dg + 2 | 3 (Dq + 12(3 

71 


+ ^^o(l - t^o) (* ■ 


(27) 



log ((Do) 


The uo/fl dependence of and D^p should be the qualitatively different for the low and 
high environment temperatures based on the spectral distribution, /((Uo) = J (cuo) coth {/3 coq) , 
of the environment. 

We have considered low (/c^T = 100 hl2), intermediate (/c^T = hl2) and high (/c^T = 
0.01 h 12) temperature limits of the environment and display the magnitude and the ratio 
of the third and the fourth order contributions to the second order one for the normal 
diffusion coefficient in Fig. 3(a)-(i). From the first column of the figure, it is obvious 

that Dxx increases with increasing initial temperature as expected, because of the increase 
in the thermal excitations with increasing temperature. Similar to the drift terms, the 
third and fourth order contributions to Dxx have opposite signs for intermediate and high T 
values while at low temperatures D^xl is negative or positive depending on ujq and 12. The 
interaction strength dependence of Dxx changes somewhat from high to low temperature 
region (compare Fig. 3(a) and Fig. 3(g); while it is almost independent of (Do = (Uo/12 ~ 1 
at high temperatures, it increases with decreasing (Dq in the low temperatures. Fig. 3(b), 

/Q\ 

(e) and (h) indicate that the ratio DxJ /Dxx is always negative and almost independent 
of the temperature. Similarly, D^]/D^] (Fig. 3(c), (f) and (i)) is positive and almost 
identical at high and intermediate temperatures while it becomes negative at low Uq and 
12 region at low temperatures. Since and D^J have opposite signs for the most of the 
considered parameter range, instead of taking into account only the second and third order 
contributions, keeping only the second order terms might be more accurate. 

To delineate the relative importance of higher order terms for the anomalous diffusion 
coefficient Dxp, we display its value and the relative third and the fourth order contributions 
as function of ujq and 12 in Fig. 5(a)-(i) at low, intermediate and high temperatures. The 
magnitude of Dxp originating from and are, also, displayed in Fig. (4) and 

Fig. (6) for different 12 values at high and low temperatures, respectively. 
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CiJq j H 

FIG. 3. Magnitude (a, d, g) and the ratio of the third (b, e, h) and the fourth order contributions 
(c, f, i) to the second order term for the normal diffusion coefficient D^x as function of the oscillator 
frequency wq and the cutoff frequency of the black-body environment If at high {ksT = 0.01 h If) 
(a, b, c), intermediate {ksT = hfl) (d, e, f) and low {ksT = 100hfl) (g, h, i) temperatures. 

Although Dxp at high and intermediate temperatures show similar behavior as the normal 
diffusion coefficient (compare the first two rows of Fig. (3) and Fig. (5), such as being positive 
for all ojq and G values and D^p and D^p having opposite signs, it becomes negative for 
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FIG. 4. Magnitude of (straight), (dot-dashed) and (dotted) contributions to normal 
{Dxx- left axis) and anomalous diffusion coefficients {Dxp'- right axis) at high temperature (/3 If = 
0.01) as function of the oscillator frequency at a) fire = 0.5, b) fire = 0.1 and c) fire = 10“^. 


a)o < 1 at low temperatures (see Fig. 5(g)). Sign of Dxp is related with the squeezing of the 
position or the momentum of the oscillator [44]. The changing sign indicates a change in 
localized dynamical variable depending on the initial temperature of the environment being 
high or low. 


The qualitative difference in both Dxx and Dxp and the contributions to them in all the 
considered orders between the high and low temperature environments is prominent when 
one compares Fig. (4) and Fig. (6). Independent of the temperature and the bath cutoff 
frequency, the shapes of the absolute values of the second, the third and the fourth order 
diffusion terms as function of the oscillator frequency are similar: the components of Dxp 
{Dxx) at high (low) temperatures display a single resonance peak at uo = fl, while the low 
temperature behavior of Dxp is similar to an anti-resonance structure. Such behavior can 
be understood by referring to the spectral density function. 
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FIG. 5. Same as figure 3 for D^p- 


IV. CONCLUSIONS 

We have considered the exact treatment of the third and the fonrth order corrections to 
the master eqnation of a charged harmonic oscillator interacting with a black body radiation 
held by using Krylov averaging method. The oscillator-bath interaction is assumed to be 
in velocity-coupling form. Diamagnetic interaction term, which is generally neglected, is 
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FIG. 6. Same as figure 4 for low temperature (/3n = 100). 

also taken into account. It is found that including higher than the second order terms 
would not change the structure of the master equation, he. they do not introduce any new 
commutators. 

For mass renormalization, decay constant, normal, and anomalous diffusion coefficients 
for the oscillator in the electromagnetic held, we obtained exact analytical expressions, to the 
fourth order, which indicate that mass renormalization and decay constant depend on the 
oscillator frequency and the bath cutoff frequency while the diffusion coefficients are sensitive 
to the bath temperature along with those two variables. Anomalous diffusion coefficient is 
found to go under a qualitative change between the high and the low temperature limits of 
the environment. 

The aim of this study was to delineate the conditions under which the higher order 
corrections to the master equation might make a difference in the dynamics of the charged 
harmonic oscillator in electromagnetic held. It is found that the single most important 
parameter in that context is the chosen cutoff frequency compared to the inverse of 
electron time Xg. We have shown that the second and the third order contributions are, 
mostly, opposite in sign and cancel each other. So, keeping only the second order terms in the 
master equation would be accurate enough for the weak interaction regime and flTe < 0 . 01 . 
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For f2 Te > 0.1, a perturbative expansion would not be appropriate as the order of consecutive 
terms have comparable magnitudes. 
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Appendix A: Krylov Method of Averaging for Deriving Master Equation 


Here, we summarize the Krylov averaging method used in Ref. [14], obtain exact analyt¬ 
ical expressions for the integrals involved in the dehnition of (ps^i) and (ps^i) aiid 
give a detailed derivation of the fourth order correction {ps,^)- 

Since the system-bath interaction is assumed to be weak, one would expect that the 
change in the product form of the initial density matrix of the total system is minimal and 
p{t) can be expressed as p(t) — ps<^ Pr where the difference from the product terms can be 
approximated as a power expansion in the perturbation parameter e and a function of the 
oscillator density matrix ps as: 

P (t) = ps it) ® pR + e (p 5 , t) + (p^, t) + ..., (Al) 

where terms containing e represent small fluctuations in the amplitude of the density matrix 
of the whole system around ps (t) ® pr. Similarly, the time rate of change of the system 
density matrix can be written as 

e iPs^ t) + ips, t) + ... . (A2) 

For the dehnition psit) = Tri:j{p(t)} to hold, Tr^lF^”^ {a, t)} = 0 for n = 1, 2,... 

. Also, B^'^^ ips,t) should be independent of bath operators, so Ttr{B^^^ ips,t)pR} = 

ips,t). 

To obtain (p^, t)s, we insert p(t) of Eq. (Al) into the von Neumann equation (Eq. (7)) 
and use Eq. (A2) for the time dependence of the reduced density matrix of the system and 
match the coefficients of the equal powers of e on the left and the right sides. The process 
produces a hierarchy of coupled equations for BP^'> ips,t) and ips,t) as: 

(9Fd) 1 

Pr (g) B^^> + = — [Hi it ), Ps ® Pr] 

Op(n) 1 

PR ® = — [Hi it ), {B^^\t) . (A3) 

m=l 

The set of equations in (A3) can be solved by starting from Rd) [ps,t) and Fd) [ps,t) and 
invoking the trace over the bath conditions mentioned above. We quote the results, up to 
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order four, from Ref. [14]: 


= (A4a) 

(fe. ii) = --R /" *2 ({12} - {1 : 2}), (A4b) 

(? h) j-oo 

{ps, h) = ^ [' dh dh {{123} - {12 : 3} - {13 : 2} - {1 : 23} 
n) J-oo J-oo 

+ {1 : 2 : 3} + {l : 3 : 2}}, (A4c) 

(ps, h) = ' dt2 / ' dts / ' dU {{1234} - {12 : 34} - {13 : 24} 

y'l fi) J— oo J— oo J—oo 

- {14 : 23} - {1 : 234} - {123 : 4} - {124 : 3} - {134 : 2} 

+ {12 : 3 : 4} + {12 : 4 : 3} + {13 : 2 : 4} + {13 : 4 : 2} 

+ {14 : 3 : 2} + {1 : 2 : 34} + {1 : 3 : 24} + {1 : 4 : 23} 

+ {1 : 24 : 3} + {1 : 34 : 2} + {14 : 2 : 3} + {1 ; 23 : 4} 

- {1 : 3 : 2 : 4} - {1 : 3 : 4 : 2} - {1 : 4 : 2 : 3} 

- {1 ; 2 : 3 : 4} - {1 : 2 : 4 : 3} - {1 : 4 : 3 : 2}} , (A4d) 

where the braces {•.. } represent the trace over the reservoir degrees of freedom of the nested 
commutators of the closed-system Hamiltonian at times ti with the density matrix ps (8) Pr 
while the column sign indicates a further trace operation: 

{1234} = Tr^{[i:f, (fi), [Hr (t^), [Hj (ts), [Hj (t4), P5 Pk]]]]} 

{12 : 34} = TTR{[Hr (ti), [Hr {h) ,TrR{[Hr (tg), [Hr (^ 4 ), Ps Pfi]]}Pfi]]}. 

Thermal expectation value of the products of held operators at different times will be rep¬ 
resented by using {ij ...) notation as 

(123 ...) = TrR{pR A (ti) A ( 12 ) A (ts)...}. 

Since A [t) is a Gaussian operator, the thermal expectation value of the product terms that 
contain an odd number of held operators are zero while those that contain an even number of 
A(ti) can be expressed as the sum of product of pair expectations by using Isserlis’ theorem: 

(1) = (123) = 0, 

(1234) = (12) (34) + (13) (24) + (14) (23). (A5) 
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Explicit form of the pair expectation (12) for the present problem can be easily calculated 
by using the time dependent held operator A{ti) as dehned in Eq. (8b) [4]: 

(12) = - duJ{u) {coth ( ^ ) cos [u (ti - t 2 )] - i sin [u {ti - t2)]}, (A6) 

TT Jo \4 Kb -l J 

where J (ca) denotes the spectral density of the bath and contains all the information related 
to the bath-oscillator interaction and is given in Eq. (6) for the present problem. 

Before evaluating B^) (p^,^ ti) for n = 1, 2, 3,4 of equations (A4a)-(A4d), we should note 
that contribution of i?A)s to the dynamics of the oscillator density matrix has the same 
form of independent of n, as: 

+ [ps{Tlp{t)+muJoT;x{t)},p{t )]), (A7) 

where Ti and T2 are time integrals of various combinations of pair correlation functions and 
can be written as Tj = 'Yhn corrections up to order n. We will provide the details of 

calculating up to n = 4. 

In the remainder of the Appendix, we will outline the evaluation of the contributions 
to the coefficients of the commutators in Eq. (A7) up to fourth order. Eirst, we note that 
(PsAi) = 0 because of the Gaussian property of A{t). The second order contribution 
is, also, easy to calculate: From Eq. (A4b), the {1 : 2} term is zero and {12} term gives 
Tf ^ = T21 and Tf ^ = T22 as: 

T21 = J *^^12 (12) COS (cjo ^12), (A8a) 

T 22 = J dfi2 (12) sin (cuo ^12) • (A8b) 

1. Third order contributions 

The contribution of {ps, ti) to the system master equation can be calculated by noting 
that only {123} term of Eq. (A4c) is nonzero. The brace {123} has a total of 64 terms; 
half of those are equal to zero because they contain odd moments. Combination of eight 
triple products of the diamagnetic term as well as eight terms formed by different 

permutations of A‘^ [ti) A { 12 ) A [t^) cancel each other. Hence, the remaining 16 terms can 
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be combined as four different integrals as = T 31 + T 33 and 


T 32 + T 34 , where 


^31 = (^ ^ 2 ^ 2 ^ dt23 Im{(12)}(23) cos (wq tis) , 

^32 = ^^12 dt23 Im{ (12)} (23) sin (wq ^13) , 

^33 = ^^12 dt23 Im{ (13) (23)} cos (cuo ti2) , 

T34 = Jo '^^ 23 lm{( 13 )( 23 )}sin(a;oti 2 ). 


(A9a) 

(A9b) 

(A9c) 

(A9d) 


One should note that when the diamagnetic term is neglected, the third order contribution 
would be zero, so the nonzero third order contributions here are due to interactions that 
contain A‘^{t 2 ) and ^^(^ 3 ). These integrals are discussed in the main body of the text. 


2. Fourth order contributions 

Evaluation of {ps,ti)j given in Eq. (A4d), is more involved and tedious compared 
to that of the lower order terms. Luckily, most of the braces in Eq. (A4d) contain a single 
element separated with columns and those are zero because of the Gaussian character of the 
bath. So, the nonzero fourth order terms can be written as: 


(psAi) 



-{14:23}}, 





dU {{1234} - {12 : 34} - (13 : 24} 


(AlO) 


where the brace {1234} has a total of 256 terms; half of those are zero because they involve 
the odd powers of the vector potential A{ti). Furthermore, 16 (48) terms which are in 
the form of product of four diamagnetic A{tiY (vector potential A(ti)^ ...) terms at times 
G, G, G; G add up to zero. The remaining terms can be written as four and six time averages 
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as: 


{1234} = ^ [(4213) [p(ti) ,p{h) psp{ti)p{t 2 )\ + (1234) [p{ti) ,p{t 2 )p{h)p{ti) ps] 

- (2134) [p(ti) ,p{h)p{U) psp{t2)] - (4321) \p{ti), psp{U)p{h)p{t2)] 

- (3124) [p(ti) ,p{t2)p{t4) psp{h)] - (4123) [p(ti) ,p{t 2 )p{h) pspiti)] 

+ (4312) [p(ti) ,p{t2) pspit4)p{h)] + (3214) [p(ti) ,p(t4)psp(t3)p(t2)]] 

+ -^ [((122334) - (221334) + (332214) - (331224)) [p (ti), p {U) ps ] 

+ ((422133) - (433221) + (433122) - (412233)) [p (ti), psP{U)] 

+ ((122344) - (221344) + (442213) - (441223)) [p (ti) ,p{h)ps] 

+ ((322144) - (443221) + (443122) - (312244)) [p (ti), p5p (tg)] 

+ ((123344) - (331244) + (443312) - (441233)) [p (ti) ,p{t2)ps] 

+ ((332144) - (213344) + (442133) - (443321)) [p (ti), psP (ta)]], (All) 

where four-time terms are due to p{t).A{t) interaction, while six-time terms originate from 
p{ti)A{ti)A^{tj).. type couplings that involve both diamagnetic and velocity interactions. 
The remaining braces in Eq. (A4d) involve only velocity-coupling terms and can be expressed 
as: 

(12 ; 34} = ^ [(12)(34) \p{ti) ,p{t2) [p{h) ,p(t4) Ps]] 

+ (12)(43) \p{ti) ,p(t2) [pspiti) ,p(t3)]] 

- (21) (34) [p(ti), [p(f3) ,p{ti) ps\p{t2)] 

- (21) (43) [p (ti), [psP {U ), p (fs)] p (ts)]], (A12) 

(13 : 24} = ^ [(13)(24) [p(ti) ,p{U) [p{t2) ,p{U) ps]] 

TTX 

+ (13)(42) [p(ti) ,p(t3) [psp{ti) ,p(t2)]] 

- (31) (24) [p (ti), [p (ta), P (t4) Ps] P (ts)] 

- (31) (24) [p (ti), [psp {U ), p (fa)] P (f3)]], (A13) 

(14 : 23} = ^ [(14)(23) [p(fi) ,p(f4) [p(f2) ,p(t3)ps]] 

+ (14)(32) [p(fi) ,p(f4) [psp{h) ,p(f 2 )]] 

- (41) (23) [p (fi), [p (fa), P (fs) Ps] P (^ 4 )] 

- (41) (32) [p (fi), [psp (ts ), p (fa)] p (f4)]] • (A14) 
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Since Eq. (All) contains both four and six index terms, one needs to transform them into 
pair product form by using Isserlis’ theorem. The formula for the four index term was given 
in Eq. (A5) while the six index term can be written as: 

(123456) = (12)(34)(56) + (12)(35)(46) + (12)(36)(45) 

+ (13) (24) (56) + (13) (25) (46) + (13) (26) (45) 

+ (14) (23) (56) + (14) (25) (36) + (14) (26) (35) 

+ (15) (23) (46) + (15) (24) (36) + (15) (26) (34) 

+ (16)(23)(45) + (16)(24)(35) + (16)(25)(34). (A15) 

In equations (A11)-(A14), we can use the time dependent momentum operator [14] 

P ih) = p{ti) cos [wo (ti - t 2 )] +muJox{ti) sin [wq (ti - t 2 )] (A16) 


to simplify the commutators involving the system momentum at different times and the 
system density matrix and obtain the contribution of 5*-^^ (p 5 ,ti) to the master equation 
Ti ^ = T 41 + Tgi + Tgs and T 2 ^ = T 42 + Te 2 + Te 4 where 


T 41 = i 
4 

T 42 = i 

Tei = 
Tq2 = 

Tqs = 

Tg4 = 




dti 2 / dt23 / dtsi {Im [(14)] (23) cos (wq tia) sin (wq t24) 


JO JO JO 

Im [(13)] (24) cos (wotu) sin (a;ot 23 ) + Ini [(14) (23)] cos (a;oti 2 ) sin (wq ^ 34 )} , 
h? w? 


I dti 2 dt23 df34 {Ini[(14)] (23)sin(a;oti3)sin(a;ot24) 

0 Jo Jo 

Im [(13)] (24) sin (wq tu) sin (wq ^ 23 ) + Im [(14) (23)] sin (wq ^ 12 ) sin (wq ^ 34 )} , 

32 a; \ roo poo 

— I) / dti 2 dt23 (if34 {Im[(12)]Im[(23)] (34)cos(a;oti4)}, 

^ ^ J Jo Jo Jo 

32a; \ 

J dti 2 J dt23 J iit34 {Im[(12)]Im[(23)] (34)sin(a;oti4)}, 

I) / dti 2 dt23 dt 34 {Im[(12)]Im[(24)(34)]cos(a;oti3) 

n m / Jo Jo Jo 

(Im [(23)] Im [(14)(34)] + Im [(13)] Im [(24)(34)]) cos (uJotu)} , 

' dti 2 I dt23 I dt 34 {Im[(12)]Im[(24)(34)]sin(a;oti3) 


32a;2 


JO JO JO 

(Im [(23)] Im [(14)(34)] + Im [(13)] Im [(24)(34)]) sin (a;oti 2 )} • 


(A17) 


Here Tgj, where i = 1, 2, 3, and 4, contain the diamagnetic interaction terms while T 41 and 
T 42 are due to p.A type interaction. Although tedious, all of T^i and Tgj integrals can be 
carried out exactly by using the rational expansion of the hyperbolic cotangent function, as 
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was the case for the second and the third order contributions. The results are: 


T41 = 






1 1 Stti 

— 7 Im [^jJl {f3, Wo)] - 7^0— csch^ {f3 uq) 


TT^ 


m 


Wo 


f-i _ _ 

I J ujn \ . r49 


' Wo 


(Al8a) 


T42 — — 
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7 ( - csch^ {(5 Wo)- 

Zi '7T 


'I3n 


TT 
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J m 


T61 = - 


Sm 

Lo^ — 
° m 


(Al8b) 
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